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FOREWORD 


Computer  programs  were  written  to  calculate  the  power  losses  of  the  dominant 
mode  on  a  dielectric  rod  at  points  where  the  rod  has  a  change  in  diameter  due  to 
either  a  sharp  step  or  a  gradual  taper.  The  equations  used  for  the  computation 
are  those  of  Dietrich  Marcuse.  The  first  program  calculates  the  power  losses  at 
a  step  in  terms  of  the  step  size  and  the  dielectric  constant.  The  second 
calculates  the  power  losses  for  tapered  roas  of  differing  lengths,  dielectric 
constants  and  taper  shapes.  The  results-  are  summarized  in  the  tables. 
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PART  A:  GENERAL  THEORY 


The  equations  governing  the  power  losses  of  a  cylindrical  dielectric  rod 
due  to  variations  in  radius  and  dielectric  constants  have  been  deduced  by 
D.  Marcuse.^  Starting  from  the  solutions  to  the  Maxwell  equations,  he 
deriv\es  the  expression  for  the  power  loss  in  terms  of  the  coefficients  of 
transmission  and  reflection  at  a  "step"  (VIZ:  a  jump  discontinuity  in  the 
radius).  These  coefricients  are  calculated  using  the  approximation  that  the 
reflected^radiation  modes  are  approximately  orthogonal  to  the  incident  modes. 

\ 

Assuming-,  cyl  inar  ical  symmetry,  the  fields  of  the  incident  bound  mode  are: 


£2  =  A  Uy  (icr)  cos  V 


Hv  =  B  J  (<r)  sin  V  4) 


Ej.  =  ”i_  ("k  a  (icr)  +  0)  JJ  B  J.^  (<r)l  cos  \) 
E^  =  i_  |So  A  ^  (icr)  +  <  0)  p  B  (Krlj 


s  in  V 


(la) 

(lb) 

( lc ) 

(ld) 


^Marcuse,  Dietrich,  "Radiation Yosses  of  the  Dominant  Mode  in  Round  Dielectric 
Waveguides,"  Bell  System  Techni\al  Journal,  Oct  1970,  p.  1665. ff. 


\  ^ 
\ 

\ 
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H,  .  -  i  [„ 


2  V 

Eq  A  _  J  ,  (<r)  +  K  3o  B  J 
2  L  °  r  °  V 


^  (< r  )J  si 


s  in  \)  ■}) 


(le) 


H,  =  -  i_  fn  K  (U  Eq  A  J'  +  8n  B  ^  J  (icr)"]  cos  V  ij)  (If) 

<j)  2  L  V  r  ^  J 


Where  \j  is  an  integer,  is  the  order  Bessel  function  and  the  factor 
exp  i(u)t  -  8^2)  has  been  omitted  from  each  term. 


k2  =  0)2  So  viq 


(2a) 


c2  =  n2  k2  -  g2 


(2b) 


n2  =  dielectric  constant 

Bq  =  propagation  constant  to  be  determined  by  the  boundary  conditions. 
The  fields  outside  Che  cylindrical  rod  are: 

E2  =  C  (iYr)  cos  v  4)  (3a) 


=  D  (iYr)  sin  V  4* 


(3b) 


^r  =  ^2 

Y-^  L 


(1)'  V  (1)  1 

Er.  =  i_  I  i  Y  3  C  H,.  (iYr)  +  0)  U  D  _  (iYr)  cos  V  $ 


(3c) 


•  r  u  (1)  (1)'  1 

ii-  S  C  .1  H  (iYr)  +  i  U)  p  Y  D  Hy  (iYr)  s  in  V  4) 
y2  L  o  r  ''  J 


(3d) 


■  r  V  (1)  „  (O’  ] 

H  =  i_  Ci)  e  C  _  H  ( iYr )  +  i  Y  3  D  H,,  ( iYr )  s  in  V  4> 

rY2L°''''  o  J 


(3e) 


H.  =  i_  i  Y  u  C  C  H  ,  '  (iYr)  +  0  D  I  H ' ’  '  ( iYr )  I  cos  V  4) 

**>  L  o 


(1)’ 

'v 


V  (1) 

’•V 


>] 


(3f) 
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=  ^2  _  1,2 

o 

To  insure  the  continuity  of  the  fields  at  the  interface,  the  following 
eigenvalue  equation  must  be  satisfied  by  3^: 


■  I  V 

<  I  j  (<3 )  <a  I  ( 1 ) 

1  /  Hi  (iya) 


-  1 


ayV-^o 


1 


tc  \  J  (ica )  <a 
1 


ya 


i  (i7a) 

o 

(1) 

Hi  (i7a) 


-  1 


(n2  -  1)  y. 

2 

< 


(5) 


(We  have  set  v  =  1,  specializing  to  the  lowest  order  mode.)  When  this 
eigenvalue  equation  is  satisfied,  the  condition  of  continuity  determines  all  the 
coefficients  in  terms  of  A.  The  coefficient.  A,  is  determined  by  the  power  in 
the  mode,  as  is  discussed  below. 


The  fields  described  so  far  are  the  bound  modes  and  they  fall  off  rapidly  with 
distance  from  rod  axis.  Inside  the  dielectric,  the  expressions  for  the  fields  of 
the  radiation  modes  are  the  same  as  for  the  bound  modes  with  A,  B  and  replaced 
by  F,  G  and  <t  respectively.  Outside  the  rod,  these  fields  are  given  by: 


9 


0  ^  H  J  (pr)  +  I  Y  (pr) 
^  V  V  _ 


+  P  OJ 

P  [k  j;  (pr) 

+  M  Y^  (Pr)j 

1  s  in  (j> 

(6d) 

Hr  =  -  i_  j 

p2  1 

0)  e  2  Fh  J 
o  r  L  ^ 

(Pr)  +  I  Y,^ 

(pr)j 

+  p  s 

(Pr)  + 

M  Y;  (pr)]j 

s  in  V  (}> 

(6e) 

"  7  1 

1 

p  U)  So  [h 

(Pr)  +  I  Y^ 

(Pr)j 

!  V 

i  +  8  - 

i 

j 

[k  (pr)  + 

M  Y^  (pr)]j 

cos  V  (j) 

(6f) 

1 

j  where  p^  = 

k2  -  02 

(7) 

The  extra  coefficients  arise  from  the  fact  that  we  can  no  longer  require  the 
fielas  to  approach  zero  so  rapidly.  Their  presence  makes  it  possible  to  insure 
the  continuity  of  the  fields  without  the  eigenvalve  equation  being  satisfied. 

1 

j  The  equations  arising  from  these  conditions  are  sufficient  to  aetermine  the 

coefficients  H,  I,  K  and  M  used  above,  and  also  to  determine  the  ratio  of  the 
coefficients  (We  specialize  to  the  case  V  =  1): 

: ^ 
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1=4/^  /  (g  -  b)^  +  (e  -  d)^  +  c~  +  V''^ 

G  j  \(g  -  b)2  +  (e  -  n-  d)^  +  +  f2  / 


where 


b  =  p/a  Jj(aa)  Yi(pa) 


(9a) 


-  1)  k  3  (£ja)  Yj  (pa) 

p  a2  a 


(9b) 


d  =  p/a  J|(aa)  J^(pa) 


(9c) 


e  =  J  (aa )  J '  (pa  ) 

1  1 


(9d) 


f  =  (n"  -Dkg  J  (ga)  J  (pa)  (9e) 

p  0^  a  ^  ^ 


g  =  J^(oa)  Yj(pa)  (9f) 

p2  =  n2  1^2  ,  g2  (10) 


This  equation  has  two  roots  corresponding  to  the  different  choices  of  sign 
for  the  ratio  F/G.  As  Dr.  Marcuse  indicates,^  this  corresponds  to  the  even 
and  odd  modes  to  be  found  in  s lab-sytnmetric  problems  (cos  <j)  and  sin  (}»  terms). 
To  completely  specify  the  problem,  we  have  to  specify  the  power  in  each  mode: 


+ 


f 


0 


(g  -  b) 


(11) 


See  footnote  1  on  page  7. 
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where  b,  c,  d,  e,  f  and  g  are  as  above. 

When  the  power  in  the  mode  is  specified;  then,  for  each  choice  of  sign  of 
F/G,  F  can  be  calculated  and  then  G  and  finally  H,  I,  K  and  M  giving  a  full 
description  of  the  incident  fields. 

When  the  incident  fields  strike  a  step,  the  continuity  of  the  tangential 
components  of  the  electric  and  magnetic  fields  requires  that 


(i)  (r)  r  (r)  (r)  1 

L  +  a  E  +  /  q  (p)  E  (p )  +  p  (p )  E  (p )  dp 
r  rr  oL’-  ^  r,2J 


c  E 
t 


(t)  r”  r  .  (t) 

+  /  q  (p)  E  (P )  +  p 

r  -b  L  t  r.l  t 


(t) 

(P)  E  (p) 
r,2 


dp 


with  similar  equations  for  E* ,  E  ,  H  ,  H.  and  H 
^  z  r  p  z 


c^  =  coefficient  of  transmission 


a  =  coefficient  of  reflection 
r 


q  (p)  =  the  amplitude  corresponding  to  the  amount  of  power  scattered  into 
the  mode,  P,  for  the  first  choice  of  sign  for  F/G 


p  (p)  =  same  as  q  (p)  with  the  sign  of  F/G  reversed. 


(12) 
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The  superscripts  i,  t  and  r  refer  to  the  incident,  transmitted  and  reflected 
waves  respectively. 

An  exact  solution  to  the  resulting  set  of  equations  would  be  extremely 
difficult.  However,  an  approximate  solution  to  the  problem  can  be  obtained  by 
ignoring  the  reflected  radiation  when  calculating  the  coefficient,  c^.  If  the 
steps  are  large,  the  radiation  will  be  predominantly  in  the  forward  direction, 
so  neglecting  the  backward  radiation  seems  justifiable  in  this  case.  If  the 
steps  are  small,  then  the  radiation  modes  of  the  reflected  radiation  are  nearly 
orthogonal  to  the  incident  radiation  modes  and  our  approximation  would  seem 
justified  in  this  case,  also. 

Witn  this  approximation,  it  is  possible  to  calculate  the  transmission 
coefficient,  c^,  and  also  the  reflection  coefficient,  ,  at  a  "step".  If  P 
represents  the  incident  power,  then: 


c 

t 


(II  +  I2)P 


r  i  1  I2 


v/here  the  two  integrals  are  given  by: 


(13) 


(14) 
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I  -  i_  (n“  OJ  €  A  -  3  B  )  (3  A  -  OJ  u  B  )  _  1  ^'^2  ^2^ 

I  ^2  o  1  11  22  21 

1  H 1  (  i  Y2  32  ' 


Ci_  +  — \j  (<  a  )  H  (i  Y  a  )  -  i_  J  (<  a  )  (i  Y  a  ) 

^2  y2  I  1  1  2  1  22  21111  21 

^*"2  ^2/  ^2 


J  (<  a  ) 

i —  (we  A  -  3  B)(6  A  -  wpB)  _ 1 _ 2 _ 2_  J  (K  a  )  H  (i  Y  a  ) 


^1  "2 


o  1  112  2  ^2  1 


1111  2  1 


H  (i  Y  a  ) 
1  2  2 


- ^  (we  3  A  A  +a)iie,B  B,) 

2  2\  0  2  2  1  *^112 


(K^  a^) 

H  Mi  Y  a  ) 
1  1  1 


J  (k.  a  )  r  I  I  2  . 

— 1 - 2 — 2 —  iy  H  (iY  a)H  (iY  a)-iY  H  (iY  a)H  (iY  a) 

1  IIII02I  20  11121 

Hi  (i  M  ®  2^  ^ 


n 

+  i_  (n^  we  3  AA+wy3  bb) 


0212  112 


2  2 

<  “  K 

1  2 


Ca.  <. 
\^2  1  1 


a  )  J  (<  a  )  -  J  (K  a  )  J  (ic  a  ) 

2  Q  2  2  o  I  2  I  2  1  i 


.  1  ■’1  ^2) 

1^2  +  y2  ^ 

1  2  Hi  (i  Y2  32) 


(a  J  (ic 
2  o  1 


1  1 
a)H  (iY  a)'-a  J  (<  a)H  (iY  a) 
2  1  2  2  1  o  1  1  1  2  1 


_LjiaJ  (k  a)H 
2  \  2  1  12  o 


(iY  a)-ia  J  (ic  a)H  (iY  a) 
2  2  1111°  2  > 
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The  subscripts  1  and  2  indicate  a  quantity  found  to  the  left  of  the  step  or 
to  the  right,  respectively.  The  quantities,  A  and  B,  are  the  amplitudes  of  the 
electric  and  magnetic  fields  inside  the  rod.  P  is  the  power  associated  with  the 
dominant  mode  (v  =  1).  With  the  aid  of  these  coefficients,  the  power  loss 
becomes : 


(17) 


This  equation  is  suitable  for  the  power  loss  computation  for  an  abrupt 
step.  Alternatively,  the  power  carried  away  in  the  radiation  modes  could  be 
calculated  directly  from: 


^  -k  /  f' 


(18) 


Both  reflected  and  transmitted  modes  are  automatically  included  by 

integrating  from  -k  to  k  instead  of  0  to  k.  If  we  take  a  =  a  (z),  then 

Aa  =  ££  Az  and  Aa  will  be  small  so  that  an  abrupt  step  is  replaced  bv  a  sequence 
oz 

of  small  steps.  The  resulting  value  of  q  is; 


q(p)  =  f  I(P,  z)  li  exp  f-i  f 
o  dz  I  J 


(6  -  8)  ds  dz 


(19) 
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The  lunccion  q  (which  governs  the  losses  of  Che  "even"  modes)  corresponds  to 
one  choice  of  sign  for  F/G,  while  p  ( p)  corresponds  to  the  other  choice  (see  the 
expression  for  F/G  above). 

The  function,  I(  p,  z),  is  calculated  by  expressing  quantities  to  the  right 
of  a  small  step  in  terms  of  quantities  on  the  left  and  a  Taylor  series 
expansion: 


F(a  )  =  F(a  ) 
2  1 


Aa  +  ... 


the  higher  order  terms  being  neglected. 


I(p  ,  z)*= _ J  (Ka)  * 

2  2  1 
4P  Y  P 


+  g)  Y  p 


A  iii  +  0)  y  B  ^ 
3a  3a 


Y  Jg  (pa)  + 
a  _ 


H  (i  Y  a) 
o 

ip  1  Ji 

H  (i  Y  a) 

1 


(pa) 


_L_  J  (pa) 
Y  P  1 


+  (6  +B)Yp(we  Aii+ojyB 
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H  (i  Y  a) 


Y  (pa)  +  1 


H  (i  Y  a) 
1 


Yi  (pa) 


y'  -p2 


-  _L_  Y  (pa) 
Y  P  1 


+  (k  +6  6) 

0 


(a  IL  +  B  m )  J 
y  3a  3a/  1 

f  "’‘■I} 


(pa) 


The  partial  derivatives  required  for  this  function  are: 


3H  =  TTP 
3a  2 


[a 

L( 


2  2 

£_  J  (aa)  Y  (pa)  + 
o  1 


L-  -  pa  + 
pa 


2  2  2 

■  J  (pa)  Y  (pa)  +  (n  P  /O  -  1)  J  (Pa)  Y  (pa)^F 
11  1  o 


2  2  1 

+  (n  -  1)  k  B  /q  j  (pa)  Y  (pa)  +  p  J  (Pa)  Y  (pa) 

U)  E  p2  p  (  o  ^  1  ° 

o 


-  i  J  (aa)  Y 
a  1 


1  '""'I 


31  ^  ^ 

Ta  2 


-  ER.  la  -  n 

2  Q  q 


£_J  (pa)  J'  (pa)  +  ll_-pa 
o  1  I  p  a 


2  2  2  J 

J  (pa)  J  (pa)  +  (n  p  /p  -  1)  J  (pa)  J  (pa)>  F 
11  lot 


(20) 


(21a) 
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+  (n^  -  1)  k~  j  (aa)  J  (pa) 


+  P  J  (aa)  J  (pa) 


0)  e  po‘ 

o 


i  J  (aa)  J  (pa)>  G 

a  i  1  I 


(21b) 


2  2 
M  =  ILS.  (n  -  1)  k 

3a  20 


(1) 


_§_L  J 

y  P  o  I  o 


(aa)  y  (pa) 
1 


+  p  J  (aa)  y  (pa)  -  1  J  (aa)  y 
1  o  a  1 


(aa)  y  (pa) 
1 


+  —1 _  J  (aa)  y  (Pa)  -  i  J  (0a)  Y  (pa)>G 

P  O  a  1  1  0  1  o  ( 


(21c) 


3M  _  _  ^  P  (n  -  1 )  k 
2a 


_ - _ <0  J  (aa)  J  (Pa) 

u)  y  p  a  )  o  1 


+  p  J  (aa)  J  (pa)  -  i  J  (aa)  J 
1  o  a  1 


^  (pa)  I  F  +|a  J 


(aa)  J’  (pa) 
1 


+  _J _  J  (aa)  J  (pa)  -  1  J  (aa)  J  (pa)>G 

P  a  a  1  1  a  1  O 


(21d) 


These  equations  are  sufficient  for  the  power  loss  computation  from  a  tapered 


rod  . 


In  evaluating  the  integral,  the  end  points  of  the  region  of  integration 

2  2  2 

where  0  =  +  k  are  singular  because  P  =  (k  -  6  )  will  be  zero  at  both  end 
points  (see  formula  18).  But  this  singularity  is  only  apparent  and  not  real  for 
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lim  I  q|  ^  =  lim  |  p|  ^  =  q 
3-*-+k  p  B-^+k  P 


The  calculation  is  tedious,  but  as  3->-k,  F/G 
modes  (F/G  >  0),  from  11,  we  can  show 


(from  8)  and,  for  the  even 


1  im  F  InP  =  p 
^-►k  ° 


exists,  while  for  the  odd  modes,  lim  F.p”^^^  exists.  For  the  even  modes: 

3->k 


-5/2  ^ 

L  =  lim  p  1(3  +  3)  Y  P  d)  e 


3h 

o  da 


J  (pa)  +  p  H*  J,  (pa)  J,  (pa) 
o _ 1 _  _  i 

.2  +2  PY  . 


Y  P 


(k  +  g  g)  J  (pa)  [a 
o  3a  1  I 


F  a  ^  a  J  (a  a)  (0  +  k)  k-Je  /y  A 

0000  0 


O  0 


a  =  Vn  -Ik  ;  F  =  lim  ^ 

°  6-^k  ^ 


-5/2  ( 

L  =  lim  p  1(6  +  6)  Y  P  d)  e  ^ 

2  (  o  o  3a 


Y  Y  (pa)  +  P  H*  Y, (pa)  Y, (Pa) 
a  o  1  _  i. 


2  +  „2 

Y  P 


Y  P 


+  (k"  +  8  g)  Y  (PanA  =  -  L 
o  3a  1  5  1 


-5/2  )  av 

L  =  lim  p  j(S+6)YPd)lJZ^ 
3  6-nc  (  o  o  9a 


Y  J  (Pa)  +  p  H*  J,  (pa)  J  (Pa) 
o  1  .1 


2  +  2 
t  P 


P  Y 


+  (k  +6  6)  J  (pa)  B 

o  1  3a 
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L 

4 


\ 


\ 


F  a(S  +  k)ka  a  J  (a  a)  5 
0  0  o  o  o 


\ 


lim  p  U  3  +  3)  Y  P  ti) 

'  o 


,,  3m  fa  Y  (P")  *  P  Yi\p 
o  3a  L  ^ 


Y^  (pa)] 

T~P  J 


2 

+  (k'  +  ft  8)  31  Y 
o  3a  1 


(pa)l  B 


L 

3 


lim  I  (p., 2)  =  Q 


For  the  odd  modes,  Che  same  result  holds,  but  its  derivation  is  easier. 
At  Che  other  end  point  of  the  region  (3  =  -k),  the  even  modes  have; 


,  .  -5/2 

lim  p  p  =  p 

3>-k  2 


while  the  odd  modes  have 


lim  F  Inp  =  p 


Then,  for  the  even  modes,  it  is  easy  to  show  that  =  0  while  the  odd 

modes  yield  the  same  result  after  much  difficulty,  although  the  calculation  is 


substantially  the  same  as  for  the  even  modes  near  3  “  • 


There  is  another  apparent  singularity  at  3  =  0  where  formula  11  gives  an 
infinite  value  of  F.  But  this.  Coo,  is  not  a  real  singularity;  for,  as  34-0 : 

CT  -*•  nk  ,  p  -►  k 
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aad  all  rhrae  limics:  iin;  F /G .  1  im  -i  / 1 3  I  F  and  1  irr  ■»  / !  3  !  G  exist.  Hence  the 

3->  o  3->.  o  ^  '' '  ' 


3-<-c 


.  ir.i  ts : 


3H 


,  lira  \^j3|  II  ,  lira  >/jiT  ||  and  1  ira  V 


3M 


’.d  1  ira  V I  3  !  all  exist  for  both. 

0  A 


choices  of  the  sign  of  F/G  and  so  Coo  does  lira  Visi  I  (P,  z).  So  the  integrand 
]3|  I~/p.  z)  has  a  finite  limit  at  3  =  0  for  both  even  and  odd  modes. 
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PART  B:  NUMERICAL  PROCEDURES 


A  computer  program  was  written  to  calculate,  from  Equation  A-17,  the  losses 
of  the  dominant  mode  at  a  step.  The  first  requirement  is  to  solve  the 
eigenvalue  Equation  (A-5)  for  Y.  The  eigenvalues  were  found  to  be  very  small  in 
the  range  of  interest  so  the  functions  of  Equation  A-5  were  expanded  in  a  power 
series  in  x  and  Inx  (x  =  Ya)  and  terms  involving  for  n  >  2  were  ignored 
(see  Appendix  A).  This  led  to  the  approximate  formula: 


Ya  ~  1.123  exp 


J 

o 


•Jl 


(<  a)' 
o 

«oa). 


(1) 


n  ' 

Where  K  =  .^n  -  Ik.  This  formula  seems  to  be  valid  as  lone  as  <  a 

o  N  ®  o 

<  2.4048,  the  first  zero  of  J  (x).  This  value  (y  )  was  used  as  an  initial 

o  e 

guess  and  a  step  of  10%  of  Y^  was  used.  The  difference  in  the  two  sides  of 
the  eigenvalue  equation  was  formed  for  several  values  of  Y  starting  with  .7  *  Y 

e 

and  increasing  in  steps  of  . 1  *  Y^  until  a  sign  change  was  produced.  A 

Regula  Falsi  method  was  then  used  to  determine  precisely  the  value  of  Y-  The 

10  2 

values  determined  for  f  =  2.75  X  10  and  n  =  2.05  are  plotted  against  k^a 
for  the  range  .7  _<  k*a  <_  3.0  in  Figure  1.  The  final  values  of  Y  are  compared 
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with  the  values  estimated  from  the  simple  formula  (1)  in  Appendix  A.  The 
agreement  is  strikingly  good  (all  things  considered)  for  values  of  k*a  <  1. 


With  the  value  of  Y  available,  B  and  K  can  be  calculated; 

o 


B/A  =  -Jls. 

(  2/ 


Boa  (i 


2  rj  (Ka)  .  1 

ica  l_Ji  (<a)  <aj 


,  i  H  (iya)  , 

_  o  - 

ya  1  ya 

.  Hi  (iya) 


This  ratio  was  calculated  and  turns  out  to  be  very  close  to  for 

.7  <  k*a  <  3.  In  the  expanded  scale  of  Figure  2,  the  differences  are  noticeable 
but  they  are  less  than  1%.  With  the  aid  of  B/A,  it  is  now  possible  to  compute  A 
from  the  formula: 


=  ^[^!!^|(aK)^  «a)  +  /  «a)]  -  2  /  (ka  )|  ^  ^ 


('■  1  \2  1 

H  (iya)\  )  2  /  u  „2\ 

_o _ j  +1  +  2S  (ka)  /l  +  B_  \ 

Hi(iya)/  ok' 


2  2  2  2  2 


\  ‘  4  /  >  J  >“0 
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Assuming  unit  power  in  the  mode,  all  the  quantities  are  known  except  A  and 

the  values  of  A  are  computed  and  plotted  in  Figure  3.  A  represents  the 

amplitude  or  the  z  component  of  the  electric  field  inside  the  red.  For  very 

small  radii,  the  amplitude  is  very  small  because  most  of  the  energy  is  outside 

the  rod  and  the  inside  fields  are  weak.  As  the  radius  increases,  the  amplitude 

of  increases  passing  through  a  peak  near  k'-a  =  1.8,  then  falling  off.  Now 

if  k-'-a  =  1.8,  then  (ratio  of  rod  diameter  to  radiation  wavelength)  will  be 

2 

about  l.S/tr  or  .57.  This  is  just  a  tad  larger  than  the  value  of  .5  which  Yip 
claims  to  be  the  optimum  condition  for  transfer  of  power  to  the  dielectric 
rod.  Yip's  work  was  for  a  dielectric  constant  of  2.56  and,  when  the  calculation 

is  repeated  for  that  value,  D/,\  turns  out  to  be  .51  (approximately). 

With  tne  aia  of  the  quantities  computed  thus  far,  it  is  possible  to  evaluate 
numerically  the  two  integrals  Marcuse  gives  and  from  zheui  to  compute  the 
transmission  ana  re  flection  coefficients  and,  thereby,  the  power  loss: 


The  values  of  these  for  various  values  of  ka  are  listed  in  Table  1. 

TABLE  1  POWER  LOSS  AT  A  STEP 


ka  j 

ct 

ar 

Ap/P 

1 .4 

.  169 

-.0207 

.971 

1 .8 

.518 

-  .0476 

.729 

2.2 

.779 

-.0649 

.389 

2.6 

.917 

-.0701 

.155 

3.0 

.975 

-.0680 

.044 

Power  Loss  at  a  Step  From  A  13,  14  and  15;  f  =  2.75  X  10 
n^  =  2.05  and  Step  Size  =  aa/aj^  =  .5. 


"Yip,  Gar  Lam,  "Launching  Efficiency  of  the  -  Surface  Wave  Mode  on  a  Dielectric 

Rod,"  IEEE  Transactions  on  Microwave  Theotj/  and  Techniques,  Vol.  MTT-18  '>‘12, 

Dec  1970. 
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FIGURE  3  AMPLITUDE  OF  THE  Z  COMPONENT  OF  THE  ELECTRIC  FIELD  (UNIT  POWER 
IN  THE  MODE)  VS  ka  INSIDE  THE  ROD 
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From  the  table,  it  is  clear  that  the  power  losses  fall  off  rapidly  with 
increasing  values  of  k'"a.  Indeed,  if  k*a  ^  1.8,  the  power  loss  varies  inversely 
with  the  square  of  k'"a  (see  Figure  4).  This  will  serve  as  an  approximate  guide 
for  calculating  losses  in  this  range.  A  complete  listing  of  the  program  is 
given  in  Appendix  B.  For  the  tapered  dielectric  rod,  the  power  loss  may  be 
computed  from  .M8-A21.  In  order  to  use  these  formulas,  the  values  of  Y,  8^,  A 
and  B  must  be  calculated  for  each  value  of  the  radius,  a.  In  the  tapered 
dielectric,  a  is  a  function  of  z  and  the  first  requirement  is  that  a  suitable 
function  be  chosen.  Four  different  types  of  function  were  used  and  the  graphs 
were  (a)  linear,  (b)  exponential,  (c)  parabolic  and  (d)  elliptic.  When  the 


functions  were  chosen,  the  axis  was  divided  by  50  points  and  the  radius  at  each 
of  the  fifty  points  calculated.  The  eigenvalue  problem  was  solved  for  each  of 
the  points  and  the  values  of  Y,  A  and  B  were  stored  in  arrays.  The  ratio  F/G 
was  evaluated  and  both  values  of  F  (one  for  each  choice  of  sign  for  F/G)  were 


calculated  and  the  corresponding  values  of  G.  The  partial  derivatives  , 

3 1  3  K  3  M 

_  ,  —  and  ~  were  then  evaluated  and,  finally,  I(p,  2).  The  ^-integral: 

da  da  da  ^ 

2 

[exp-i  /  (6^  -  6)  ds]  was  evaluated  using  Simpson's  rule  when  z/Az  was 

even.  When  z/ Az  was  odd,  Simpson's  rule  was  used  to  calculate  the  integral  up 

to  z-Az  and  the  trapezoidal  rule  to  extend  the  integral  to  z.  With  I(p,  z)  and 

the  8-integral  both  defined,  Simpson's  rule  is  used  to  calculate  q(p)  &  p(p) 
from  A19  and  used  again  to  calculate  the  power  loss  from  A18. 


For  the  linear  tapers,  the  step  was  always  .5  and  the  length  was  varied  to 
make  the  taper  more  or  less  steep.  The  dielectric  constant  was  also  varied  to 
study  Its  effects  on  the  losses.  The  results  are  summarized  in  Table  2. 
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TABLE  2  POWER  LOSS  OF  A  TAPERED  ROD  (LINEAR  TAPER) 


* 

1.5 

2.05 

2.56 

t 

10 

.42482 

.11898 

.034436 

L/a  ^ 

20 

.30241 

.067682 

.016090 

100 

.045451 

.010152 

.002784 

Power  Loss  (AP/p)  as  a  Function  of  Dielectric  Constant  (n^) 
and  Length  to  Radius  Ratio  (L/c  )  for  f  =  3  X  10^^, 
a^  =  .4  X  10~^,  a^  =  .2  X  ^a  ^  =  2.513. 


As  is  readily  seen,  long  rods  lose  less  power  than  short  rods  and  high 
dielectric  materials  lose  less  than  low  ones.  For  exponential  tapers  with 
exponent  4.6  (ie:  exp  -  4.6*  z/L),  the  results  are  summarized  in  Table  3. 


TABLE  3  POWER  LOSS  OF  A  TAPERED  ROD  (EXPONENTIAL  TAPER) 


•n2 

1.5 

2.05 

2.56 

10 

.45190 

.12191 

.020800 

L/a^ 

20 

.36540 

.045415 

.0025997 

100 

.054494 

.0012596 

.0022469 

50 

.0016464 

Power  Loss  as  a  Function  of  Dielectric  Constant  and  Length 
to  Radius  Ratio  for  an  Exponentially  Tapered  Rod. 
(Frequency  &  step  size  as  in  Table  2) 


Comparison  of  Tables  2  and  3  show  that  the  exponential  taper  is  not  verv 
much  better  than  the  linear  taper  in  the  region  considered.  The  question 
naturally  arises  as  to  whether  the  exponent  -4.6  is  the  best  choice.  So  a 
family  of  exponential  curves  was  chosen  (exp  ~  p  *  z/L)  and  the  dielectric 
constant  was  taken  to  be  2.05.  The  calculations  were  then  performed  for  various 
values  of  p.  If  p  -*•  «,  the  curve  approaches  a  step  and  the  results  of  the 
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preceeding  section  are  applicable.  If  p  -►  0,  the  curve  approaches  a  straight 
line  and  the  linear  taper  calculations  apply.  The  results  are  sununarized  in 
Table  4. 


TABLE  4  POWER  LOSS  OF  A  TAPERED  ROD  (VARIOUS  EXPONENTS) 


L/a^- 

3 

10 

30 

100 

taper 

0 

.119 

.010152 

Linear 

3 

.227 

.110 

.02437 

.00131 

e~3z/L 

+ 

p 

4.6 

.122 

.00126 

g-4.6z/L 

4- 

7 

.244 

.163 

.02364 

.00230 

e-7z/L 

00 

.194 

.194 

.194 

.194 

STEP 

Power  Loss  for  Various  Exponential  Curves  and  Various  Lengths. 

n2  =  2.05. 


From  the  table,  it  appears  that  minimum  power  losses  occur  for  different 

j  exponents  depending  on  the  value  of  L/a^-  ^or  L/o^  =  10,  the  minimizing  value  of 

i 

'  p  should  be  close  to  3,  while  for  L/a,  =  100,  it  should  be  nearer  to  4.6.  The 

I 

I  dependence  of  p  on  the  dielectric  constant  would  require  further  calculation. 

! 

Two  more  tapers  were  calculated  for  various  values  of  L/A.  The  first  was  a 

parabolic  taper  with  the  vertex  at  the  thin  end  of  the  rod  and  the  second  was  an 

I  ellipse  whose  center  was  above  the  thin  end  of  the  rod  and  to  the  left  of  the 

1 

i  thick  end  (see  Table  5). 

I 
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TABLE  5  POWER  LOSS  OF  VARIOUS  TAPERS 


Taper 

3 

10 

30 

100 

Ellipse 

.21935 

.10793 

.021004 

.0015239 

Parabola 

.229 

.11975 

.036586 

.001378 

The  elliptical  taper  has  lower  losses  chan  the  parabolic  for  short  antennas 
and  Che  same  is  true  when  the  elliptical  rod  is  compared  with  the 
exponential  and  linear  tapers. 

A  complete  listing  of  this  program  is  given  in  Appendix  B. 
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APPENDIX  A 

SOLUTION  OF  THE  EIGENVALUE  PROBLEM 


For  the  bour.d  nodes,  it  is  necessary  that  3^  satisfy  the  eigenvalue 
equa  t ion  . 


i  (i'l'a) 
0 


2  /J  (<a) 


Hi  (iY 


_ -  n  £L  /  _£ _  -  i_ 

<  I  J  (<  a )  <  a 

iYa)  'I  '' 


Ya 


i  (iYa)  ^2  /J  «a)  , 

o'  _  aY  /  o  _  1 


Hi  (iYa) 


I  J  (<  a  )  <  a 

\  1  > 


/  2  2  1 

1  2  2  2 

J  2  2  2 

wnere  Y  -  ' 

^'3  -  k  and  <  = 

n  k  •  3  = 

0 

\(n  -  1)  k  -  Y 

Set  ya  -  ^  >  Ka  =  iC"  ana  ka  =  k* 

J  (k*) 

Expand  the  term  — _ _ -  -  __  in  a  power  series  about  x  =  0  and  set  ti 

J ,  (< " )  ,  rr - 

t  =  *-  ■  a_  .  where  ^*  =  '\n‘-  ~  I  k* 

1  J  (<",'  <■'  o 

1  o  o 


A-L 
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(= 


+  2  In  +  (n  +  1)  I  = 
7 


2  4 

_2 _  +  0(x  ) 


2  ••»2 
o''  k  (n  -  1)  k 


+21ni+(n  +1) 


2  Y  +  2  In  2.  + 
1  2 


(n^  +  1)  t. 


/  2 

2  In  2i  = 

II 

I 

r 

? 

+  1 

A 

2 1  2 

k  yn 

-  i 

II 

1 

r  1  n^  + 

i  + 

*2  2 

k  n  -  1 


2  Y. 


A.  . 

1  +  _ 

2 

.1 

2 

J  k 

L  (n  -  1 )  J 

_  \_ 

r„j  *  ii 

0 

: 

^  1 

k 

L  n 

-  1- 

0 

Z. 

-  n 

+ 

L  t 

k'"'’ 

1 

0 

+  1 

j 

0 

(K-) 

7 

_  n*"  +  1 

ii 

< 

0 

J 

1 

(<*) 

0 

*  2 

K  ^ 

0 

2 

+  0(x  ) 


*  2  ^ 
K  =  (n  ^  1 )  k 
0 


7  2  1  I 

n*'+l=n  +1  1_  =>ln^=“_<2Y  + 

^,v  2  2  ^-2  2  2  I  1 

o  n  -  1  k 


2  ,  J  (k*)) 

n  +  1  o  o  ( 

a  *  I 

ic  Ji  (iCq)  / 

rt  ^ 


In  21 
2 


2  ,  1  J  <<*) 

Y+n+11  o  o 
1  2  * 


o  Jl 


+  0(x  ) 


X  a  2  exp  (-Y  )  exp 
1 


(2  ,  1  J  (<*)  ) 

^  1  1  o  o  I 

2  V*  *  ( 

o  Ji  (<o)^ 


A-3 


SSWC  TR  81-299 


,  2  ,  J  (iC  a) 

Ya  -1.123  exp  -  i  P  -2 - 2__ 

,  2  <  a  J  (k  a) 

o  1  o 


This  formula  gives  an  estimate  for  the  value  of  y  that  satisfies  the 

eigenvalue  equation.  The  estimate  will  be  reliable  when  (1)  Ya«l  (2) 

I  2  I 

V  n  -  1  ka  <  the  first  zero  of  J  =  2.4 


TABLE  A-1  SOLUTION  OF  THE  EIGENVALUE  PROBLEM 


For  n2 

=  2.05 

ka 

<2.4  =  2  34 

^  1755 

ka 

Ya  (formula) 

Ya  (computer  calculation) 

• 

.5 

1.5  X  10-5 

1.48  X 

10-5 

.625 

9.74  X 

10-^ 

9.76  X 

10-5 

.75 

9.45  X 

10-3 

9.50  X 

10-3 

.875 

3.77  X 

10-2 

3.77  X 

10-2 

1.0 

1.08  X 

10-1 

9.32  X 

10-2 

1.125 

1.69  X 

10-1 

1.25 

2.63  X 

10-1 

1.375 

3.56  X 

10-1 

1.5 

4.72  X 

10-1 

5,48  X 

10-1 

For  the  range  of  values  in  Table  A-1,  the  formula  seems  to  be  accurate  to 
within  15%.  It  gives  a  useful  initial  guess  for  the  start  of  an  iteration 
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APPENDIX  B 

FORTRAN  LISTING  OF  THE  COMPUTER  PROGRAM  USED 
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PROGRAM  TAPER(IKPUT,0UTPUT, TAPc6> 


common 

A. 

BETA, 

3ETA0, 

BINT, 

CAPA, 

CAPS, 

2ChCA(Z),CI0A(2) 

,CKnA(2l  , 

DMOA  (2)  , 

EPSO  , 

Ft?)  , 

FOGt?) 

TFUNC, 

G(?l , 

gamma. 

IfiZl, 

IRI2, 

JOR, 

JOS, 

4  JIR, 

JIBP, 

JtS, 

JISP, 

KAPPA, 

KSC, 

SLENGTm, 

MU, 

HU1, 

NSC, 

OMEGA, 

Pf 

PI, 

e  rhq. 

sigma, 

YOO, 

YOS, 

Y1R| 

YIRP, 

YIS, 

7  YtSP 

cc“MON  /number/  rma)(,nmax,ift,ipthi,ae,alen, switch, Z,0Z 

COMMON  /*fiY/  AARRAV 

ccmmcn  /cry/  BARRAY 
COMMON  /GRY/  GARRAV 
CCmmon/BIM/BInTI 
dimension  eiNTKSl) 

DIMENSION  GAPRAYC  51  ),AARRAYI  51  )fBARRAY(  51  ) 

real  LENGTH, N,KSC,NSO,LOA  »MU,MU0 ,NU, <APDA 
C  all  dimensions  must  be  MkS  StStlSSSflSSSSSSttfttttftSSSSSSStS* 
p»l. 

PIm3. 14155 

"Uil*4.E-7*PI  S  EPS0=l.E-5/36./PI 
IPT  314,4  t  IPTMl  =IPT-1 

N“AX  *  50  I  KMAX  =50 

KMAXl  =  NM4X  *1 

3  READ  5,A,.IOA,STEP,NU,NSO  , SNITCH 

5  FORMAT (6E12. 51 
IFCLOA.LE.O.)  STOP 
OMEGA  =2.»P1*MJ 
AE  »  STEP*A 
LENGTH  =  ICA  ‘A 
PRINT  6,A,AE,LOA,NU,NSO 

5  FCRMAT(*1THE  GFCMFTRy  !$•,/,•  EEGINNING  RADIUS  =  •,£!?. 5, 

C  •FINAL  RADIUS  =  •,El?.5,/,«  LENGTH  =  •,E12.5,//,»  FREQLENCY  =  *, 

A  E12.5,»  DIELECTRIC  CONSTANT  =  *,E12.5) 

ALEN  =  length 
MU»  INITflJ 
MU3MU0 

KSO  =OMEGA»OMEGA»MUO»EPSO  I  K=SORT(KSO) 

CALL  GAMSGLV 
CALL  eiNTC-RL 
DINT=0. 

SIMPSON  =4.  t  SIGN  =-l. 

beta  =-K  t  08ETA  =K/FLOAT{KNAX  ) 

KMAX  =  2*  UMAX  -1 
DO  100  KKxl,  KMAX 

BETA  =  BETA  ♦  DEETA  R  RHO  sSORT CkSO-BETA*BETA) 
IF(ABSIBETA).GT.0.05<OBETA)  GOTO  90 

CALL  6ETEQ0 
DINT  =  DINT*  FUNC 
GO  TO  99 

90  SIGMA  =  SCRT(NSQ»KSQ-9ETA-8ETAI 
CALL  INTEG 

DINT  =  DINT*  SIMPSON»FUNC*AeS(BETA»/RHO 
99  SIMPSON  =  SIMPSON  ♦2.»SIGN 
100  SIGN  a  -SIGN 

DROP  =  DINT*CBETA/3. 

PRINT  l2,OINT,DeETA,OPOP 

12  FCPMATCIHI,  E12.9,*0BETA=»,E15.9,*AN0  THE  PONER  LOSS  IS», E 12. 5, // ) 


KMAX  = 

tKMAX  ♦U/7  S  GOTO  3 

STOP  S  END 

SU9R0UTINE  BINTGRL 

COMMON  A,  beta. 

SETAO, 

BINT, 

CAPA, 

CAPE, 

20M04  r?) 

lOIDA  (?)  ,0K0A(2)  , 

OMOAtZ)  , 

EPSO  , 

Ft?)  , 

FOGt?) , 

3FUNC, 

Gt2l,  gamma. 

IRZl, 

IPZ2, 

JOR, 

JOS, 

4  JIR, 

JIRP,  JlS, 

JISP, 

KAPPA, 

KSC, 

Copy  available  to  DTIC  does  not 
pezxnit  fully  legible  tepioductioii 
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SlENGTH,  HU, 

MUO, 

NSQ, 

OMEGA, 

P, 

PI, 

6  RHO,  SIGMA, 

YOR, 

VOS, 

YIR, 

YIRP, 

YIS, 

7  YISP 

COMMON  /NUMBER/  KMA *. NHAX ,IPT , IPTMl , AE , «LEN, SHI TCK, 2,0 Z 

COHHON  /bin/  aiNTl 

CCmhCN/GRY/GARRAY 

DIHENSION  GAPRAYIl  I.BINTHl  I 

«EAL  IRZlf IRZ2,JOR,JOS,J19, J1S,J1PP, JISP, K, KAPPA , KSQ, HU ,HUO ,NSQ 
peal  LENGYH 

USES  A  CCH9INATI0N  OF  SIMPSON'S  PULE  AMD  THE  TPAPEZOIOAL  R6LE  TC  EVALU 

THE  INTEGRALS 

BINTHDaO. 

FIRST  =  SQRKKSQ'GARRAYCU'GARRAYril  ) 

DC  5  MNaZ.NHAX  ,2 

second  =  ?QRT(KSQ*  GARPAY (NN) 'GARRAY tNN)  ) 

BINT1(NN»  =  BINTKNN-U  ♦  .5«(FIRST'SEC0N0) 'DZ 
NNP  =  NN*1 

SECOND  *  FIRST  *4 .'SECOND 

FIRST  S  SORT(KSQ'  GARPAY <NNF» 'GARRAY (NNP)  ) 

5  aiNTirNNPI=  BINTKNN-l)  ♦  CFIPST'SECONO  > 'DZ/3. 

RETURN  t  END 
SUOROUTIME  INTEG 


common  a. 

beta. 

betao. 

BINT, 

CAPA, 

CAPS, 

20H0A(2),CICA(2) 

,OKOA(2»  , 

0M0A(2), 

EPSO  , 

F«2)  , 

F0G(2) , 

3FUNC,  G(2J, 

gamma. 

IfiZl, 

IRZZ, 

JOR, 

JOS, 

4  JIR,  JIRP, 

JIS, 

JISP, 

X, 

KARRA, 

KSC, 

slength,  mu. 

MUD, 

NSQ, 

OMEGA, 

P, 

PI, 

6  RHO,  SIGMA, 

YOR, 

yOS, 

YlR, 

riRp, 

YIS, 

7  YISP 

COHNON  /VUHBEfi/  KHAX,NHAX,IPT, IPTM1,AE, ALEN, SWITCH, Z,OZ 
CCHHON  /BIN/  BINTt 
DIMENSION  BINTKU 

REAL  JOR,JQS, Jlfi,J19o,JlS, J ISP, K, KAPPA, KSO, LENGTH, MU, HUOjNSQ 
REAL  IfiZlilRZ? 

CChplEX  A»G,SUN1,SUH? 

PRINT  2,  beta 

2  FCPHATC  IHl,'  BETA  «',  E12.5,'  SIttSStStttSSS*  > 


IPOT=l 

NHAX  1=  NPAX  *1 
Z-<! 

A*FN(7I  t  CAOZ*  FNPtZJ  J  P=l. 

CALL  GAHFINO 
CALL  ERROR 
CALL  IRZC 

PRINT  3,A,Z,IfiZl,IRZ’ 

SU«1  *  IRZt'OAOZ'BEJl (KAPPA'AI  /GAH"A/GAMMA 
SUH2  »  IRJZ'OAOZ'PE  JUKAPPA'A) /GAHMA/GANMA 
SIMPSON  »A.  1  SIGN  =-l. 

OZ  a  LENGTH/FLOAT (MMAXI  t  Z=  0Z7NHaXMl=NHAX-l 
DC  100  NN=1,NMAXH1  SA=FN(Z)  t  OAOZ=FNP(Z) 

CALL  GAMFINO 

BINT  a  BINTKNN'll  -RETA'Z 

ARG  a  CHPLX<0.,EINT1  *  ARG=CEXPC-ARG) 

CALL  IRZC 
CALL  ERROR 
IF(IPPT.LE.l) 

IPRINT  3,A,Z,IRZ1,IRZ-',6INT  , GAMMA  ,CAPA,CAPe 
3  FORMATdH  ,'A  a' ,  E 1 2.5  , '  Za' ,F1Z.5,'  '  ,  6E12.5) 

IFdPRT.GT  .IP  TPl)  lORT  a  IPRT  -IPT 
IPRT  a  IPRT  ♦! 

SUHlaSUnlA  IRZl'OAnZ'ARG'SlHPSCN  'PEJT ( KAPPA'A > /GAMMA/GAMMA 
SUM2aSUM2*  IRZZ'OAOZ'ARG'SIMPSCN  'BEJl (KAPPA'A)/GAMMA/GAMmA 
SIMPSON  a  SIMPSON  'P.'STGN  t  SIGN  =-SIGN 


a 

o 

A  9 

S 

r\  O 
U 


<a 


ing  7sY»0Z 
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A=CH{Z»  J  0*0Z  =  FNP(ZI 
CAL'-  GAMFTNO 

BINT  =  BINTltNMAXH  -  BFTA*? 

ARG=CMPLX«0.,elKT)  »  APG  =CE)'Pt-AR6) 

CALL  IR2C 

PRINT  3|A,Z,IR2l,IRZ? 

SUMi=  SUP1»IOZ1»OAOZ»*RG  •BEJ1«'<APPA»A)/GAMMA/GAHMA 
SUM2=  SUH2*IP??«0A0Z»ARG»9EJ1(KAPPA«AI  /GAMHA/GAHNA 
FACTORS. 26’PI  /RHC/RHO 

SUN1=  SUP1»FACTCR  *  SON2=SUH2*FACTOR 
SUNlsSUPl»0Z/3.  S  SUM?3SUM?»0Z/3. 

FUNC=  CAB?  (SUMl) ’‘Z  ♦  CABS(  SUNZ)  *•  2 
PRINT  4,FUNC 

4  FCRHATC  THE  VALUE  OF  0  SOUAREO  ♦  P  SQUARED  FOR  THIS  VALUE  OF  BET 
TA  IS»,  £20.12) 

RETURN  t  END 
SOROUTINl  IRZC 


common  a. 

beta. 

BETAS, 

bint, 

CAPA, 

capb. 

20HOA(2),010A(2) 

,CK0A(2) , 

DMOA  (2)  , 

EPsn  , 

Pt2»  , 

F0G(  2) 

3FUMC,  C-(2), 

gamma. 

IRZl, 

IRZ2, 

JOR, 

JOS, 

4  JIR,  JlRP, 

JIS. 

JIS®, 

KAPPA, 

KSC, 

5LENGTH,  MU, 

MUO, 

NSQ, 

OMEGA, 

P, 

PI, 

6  RhO,  SIGMA, 

YOR, 

YDS, 

YlR, 

YIRP, 

VIS, 

7  VISP 

OIhEnSION  H(3),  HF(3) 

OINENSION  IRZ(2) 
real  IRZ1,IR22 
real  IRZ  ,mu 

real  JOS.JlSi JlSPfJOR .J1R»J1RP*NSQ, k,muo, mu, KSQ, kappa 
call  HANKEL<1,GAHMA*A,MF,H,F0,H) 

CALL  FOGC 

TEOM  =  A»«GAMMA»JOR-phC*EO/HF(1).jip)/(GAMMA»GAMMA*RHO»RHC) 

TEOHs  TERP-JIR/RHC  /GAMMA 

TEOmi*  A*(GAMMA»YOR-RhO»F9»T1R/MF(1)  ) 

TE»M1=TE RMl/ (GAMMA* gamma* OHO»RHO) -TIR/GAMMA/RHO 
CALL  OHOAC  I  CALC  OkOAC 

call  oioac  s  call  QMOAC 

FACl=(8ETA0*eETA) ‘G AMMA*RmO*OM£GA 
FAC2*  FACI*MU*CAop*FACl=FACl*EFSO*CAPA 
FACl=  FACl*TERM*(KSO*BETAf)*SETA)  •CAPB'JIR 
FAC2=  FAC2*TERM*(KSO*9ETAO*BETA)*CAPA  'JIR 
FAC3=(eETAO*BETA)*GAMMA*RHO*OMEGA 

FAC4=FAC3*MU*CAP3*TERM1  S  FACZ*  FAC3*EPS0*CAPA*TERM1 
FACTS  FACI*rifi*CAP8*(KSO*9ETAO*0ETA) 

FAC4=  FAC4*Y1R*CAPA*(KSO»BETAO*9ETA) 

00  2  11=1,2 

2  IR2(III  =  FAC1*CHCA  (II)  *FAC2*CK0At  II)  ♦FAC3*0IDA(  1 1 )  ♦£  AC4  *CMO  A  ( 1 1) 
IR2I  =  IRZd)  1  IR22=  IRZ(2) 

RETURN  *  END 
SL'3RnUTINE  BETEOn 


COMMON 

A, 

beta. 

3ETA0, 

bint, 

CAPA, 

CAPB, 

2nHnA(2),0I0A(2) 

,OKOA(2)  , 

OMOA  (2)  , 

EPSO  , 

F(2)  , 

FOG(  2) 

3FUNC, 

G(2)  , 

GAMMA, 

IRZl, 

IR22, 

JOR, 

JOS, 

4  JIR, 

JtPP, 

JIS, 

JISP, 

K, 

KAPPA, 

KSC, 

5LENGTH, 

MU, 

MUO, 

NSC, 

omega. 

P, 

PI, 

6  RHO, 

SIGMA, 

YOR, 

YOS, 

YlR, 

YIRP, 

VIS, 

7  TlSP 

CCMMCN  /bin/  BINTI 
dimension  SINTl(l) 

COMMON  /NUMBER/  KMAX,NMAX,IPT,IPTlil,AE,ALEN,SHITCH,Z,OZ 
"EAC  length, nSQ 

real  IRZl,  IRZ2,J0R,J'»S,  J1P,J1S,  J1RP,J1SP,K,  kappa,  KSQ,  mu,  MUO,  NSQ 
CCMpLEX  ABG,SLMl,SUM? 

PRINT  2,  IETA 


.  DTlC  does  oo' 
Qvoilable  to  D _ 
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7  FC°“ATtlHi,»  EETa  =  »,El2.S,*  !  JJ T *  5  i  «!$$??$ r* ) 

IFOT=l 


?£T4  =f). 

STC-Mfl  =K»5CRT(NE0) 

2  =0. 

A  =  rN<Z)  t  0402=  F^JPJZ* 

NMAXI  =  RHAX  ♦! 

C4LL  G4»<F1\D 
CAUL  IRZO 

PRINT  3,  «i  2,1=21,  IPZ’ 

3  FCOHATdH  ,•4=  •,F12.5,»Z=  •,E12.5,* 

sum  =19  21  •DAC2*eEJl(<APp4»  A)/;A«‘*4/GAN,,'4 
SU'“2  =IR  22  •0  4C2*9E Jit  KA=P 4*4)  /GANPA/ GAM-4 
SI-PSON  =9.  S  SIGN  =-l 
2  =  02 


NH4X-1  =  h-4X-l 
CO  100  NN  =1 ,NP4XM1 
4  =  FN(2)  t  0402  =  FnotZ) 
call  GAMFInO 

3INT  =  eiMl  (NN»  II 

4RG  =CMPLX(0.  ,SINTI  t  49G  =  CEXP(-4RG) 

C4LL  1920 
C4LL  ERROR 
IFtIPPT.LE  .1) 

loeiNT  3,4,2,IR2l,  IR2->,pINt,  G4«“4  ,  C  4  =  4  ,  C  4P- 
I  F  ( I  PR  T.  GT  .  I  P  Tfll  IP9T  =  IPRT  -IPT 
IFOT  =  IPfT  *1 

SU'«1  =  SU-l  ♦  IP  21*  Cl  02*4pG*SIPPSCN»  =  EJ1  (<4PPt  *4)  /G4M-4/G4HM4 
SU-?  =  SUM?  ♦  I  R  Z7»  C102*i  PC- ‘SI-PSCN'SEJl  t  ><499  4  I /GiM-4/G4HM4 


100 


SI-PSCN  =  SIMFSON  ♦?.»SIGN 
SIGN  =  -SIGN 
2=  2  ♦  02 

4  =  FNtZI  I  0402  =  FNP{2) 

C4LL  G4NFIN0 

SINT  =  eiNTKNMAXll 

4RG  =CHPL>(0.  ,eiNT)  I  4RG  =  CEXPt-4RG) 

C4LL  1920 

PRINT  3,4,2,TR21,IP2-> 

SL'“1  =  SUI*1  ♦  IP2l*C402*4PG  •BEJlt<4PP4*4)/G4M“4/G4MM4 

SU''7  =  SU'f?  ♦  1  =  27*0402*490  •EE  J  1  t  <  4  PP4  •  4  I /G  4M- 4 /GAH-fl 

r4CT0P  =  .  7G*FI/RI'0/PWC*02/‘'. 

SC-l  =  SU'NdFACTQP  *  SUM2=  SUP2*F4CTCR 
FONC  =  CAiStSUMl)  ••?  ♦  CAPS  ISU’'?)  *•? 

FUNC  =  FUNC/P.HO 
return  J  ENO 
SU’PROUTINE  BESSFN 


CC--CN 

/G4 “SOL/NMAX 

CC-MCN 

A, 

seta. 

3ET40, 

BINT, 

CAPA, 

OAFS, 

70HO4 (7) 

,0104  (7) 

,0^04(71  , 

ONOA  (71  , 

EPSO  , 

F(2)  , 

F0G(  71  , 

IFUNC, 

G(2)  , 

GAM- 4, 

I  RZ 1 , 

IP2?, 

JOR, 

JOS, 

R  JIR, 

J1PP, 

Jts, 

JlSP, 

K, 

XAPOA, 

ESC, 

5LENGTH, 

-U, 

MUO, 

NSC, 

omega, 

P, 

PI, 

6  9-0, 

7  YISP 

SIG-4, 

Y09, 

YCS, 

YIR, 

YIRP, 

YIS, 

9E4L  rR21,I92?,JOP,JOS,Jl=,JlPP,JlS,JlSF,<,KAPPA,<SO,LENGT-,NU,PUC 
9ESL  NSO 

o  =  Rh0*4  «  E=  SIGNS*® 

J0==?EJ0<R)  f  T09==ET0(P> 

JlP=9FJlt  =  )  T  rt==E'‘Tl(R) 

J1=P  =  IR*J0R-J1P)/P  S  T1=P  =(R*T0R-TlPl/P 


J0S=BRJ0(SI  *  YOS=3FrO(S) 

JlS=flEJl(Sl  I  >1F=°FV1(S) 

JISP  =  ( S* J0S-J1S)/S  *  NISP  =($*Y0S-viS)/S 
RETURN  I  ENO 
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cUSRQ'JTINf  IPZO 

nOMHON  a,  eET4,  3£TaD,  EIM,  CAPi,  CAoa, 

3nv-’'4(2),CIEA(?)  jOKOSI’*  ,  0^04(2),  CPSQ  ,  F(?),  rOG(Z), 

3FUNC>  G  131  f  GAKM4,  IPZl,  IRZZ,  JOR.  JOS, 

A  JIR,  J19P,  JlS,  JISP,  K,  kappa,  KSC, 

5LSNGTH,  MU,  MUO,  NSC,  OMEGA,  9,  PI, 

6  RMO,  SIGMA,  V09,  YOS,  YIR,  Y19P,  Yl$, 

7  YISP 

COMMON  /NUMBFR/  KMAX, NHAX ,IPT, IPTM1,AE, ALEN, SWITCM, Z,0Z 
dimension  mF<3),CUM(7)  ,  IRZ(?) 

real  JOP , JOS, Jlfi, JlS, J1°P, JISP, IRZ,K, KSO,HUO, NSC 
PEAL  IRZl, IRZ?, kappa, length, MU 

PHO  =K 

SIGMA  =K»SCRT (NSCl 

CALL  eSSSPN 

a  =  K»Yir» JISF/SIGMA 

0=K»J1R«J1SP/SIGMA 

E=J1S*J1RF 

GG=JlS»YlRo 

Fl=  SORT  (PUO/EPSO) 

FCG(l)  =  <GG-BI»  (GG-3)  ♦(E-0)*(F-D) 

FCG(l)  =F0G(1)/ ( (GG-MSC*B)»*?  ♦(E*NSC»D)**7) 

FCG(l)  =SQRT (FOG ( 1)  1  ‘FI 
FCG(2)  =  -  FOG(l) 

FAC  =  .125*PI»»3»A»A»OHEGA»£PSO/RhO 
T1  =  (GG-NSQ»£1  *•?  «■  (E-NSQ»0)»*Z 
TZ=<GG-a)»‘?  ♦  tE-0)»*2 
T?  =  Ta-'MUO/EFSO/FOGfll/FOGU) 

Ftll  =  FAC*(T1  YT?! 

F(t)  =S0OT  (l./F(l)  ) 

F(?)  =  F(|) 

G(l)  =  F (1)/FCG(1) 

•G(?l  =  -  Gtl) 

OMOA(l)  a  PI»K».?»(NSO-l.l*'<»A»Ptl) ‘JlS 
CIOAllI  =  -0HCA(1)»J1R 
OHPAd)  =  OHOAdloYto 
OmOA(?)  =  CHOAd) 

OIOAC?)  =  OIDAd) 

DK9A  d)  =.5»PI* (NSa-l.l*KSa/SQRT (NSO) 

Tl  =  A»  JOS 

TZ  =  ?.•  J1S/A/KSG/SCPHNS0) 

T3  =  -JlS/K/SCPT  <NSn) 

DMDA(l)  =  -DKCA  dl*  (Tl»JlPP+T2»JlRfT3*JCR)»Gd) 

0M94(2I  =-OMOAd» 

OKOA  d 1  =  CKCA ( 1) * (Tl»YlRP+TZ»r IR+TZ’YOS) ‘G  d) 

OkOA(?)=  -OKOAd) 

Tl  =9£rAO*GAMMA»K»OM£GA 
TZ  =  Tl»MLO»CAPe 
Tl  =  T1»EPS0*CAPA 
G4MMA»A 

CAUL  H4NKELd,X,HF,DUM,F0,O'JM) 

T£0M  =  gamma  •JOR -RHC*FO/HPd)  ‘JIR 
T£PM  =A«T£RM/ (GAMMA»GAHMA  ♦  RHC'RHOl 

Term  =  tefm-jir/rmo/gamma 

Tl  =  T1»T£RM  S  T2  =  TZ’TERM 

Tl  =  Tl  ♦  KSC*C40p»J1P 

TZ  =  TZ  ♦kSQ»CAFA»JIO 

TFPM  =GAMM  A»y  0R-PmO»F0/HF  (1  )»YiP 

TEPM  =  T£R“*A/(GAHHa»GAMMAYKSC)  -YIP/GAMMA/K 
T3  =PETA(J»GAMMA»y»0“£SA 

TA  =  T3‘MLa»CAPE»TCPM 
T3  =T,3»EFS0*CAP  1»T£PM 
TZ  =  T3  fKSQ’YlR’CAF? 

Ta  =  Ta  ♦KSC»riF ‘CAPA 
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DC  10  1=1,? 

10  =  (Tl»0HD4tI)  ♦T?»DK3a(T)  ♦  T  3»  DIDA  ( I  )  ♦!<,  *OMD«  (  1 1  ) 

IPTl  =  IP2(1)  ♦  IPZ?  =  IR7l?> 

RETU9N  i  END 
*:U9ROUTINF  FQCC 


CCNMON 

A, 

BETA  , 

BETAO, 

BINT, 

CAPA, 

CAP?, 

20HnA (?) 

,D10A(?) 

,DKOAC>)  , 

DMCA  (?)  , 

EPSO  , 

H(  2)  , 

FCG(’ 

3FUNC, 

Q  (21  , 

G  A  MM  A  , 

IRZl, 

IRZ?, 

JOR, 

JOS, 

>*  JIR, 

J  IRP, 

JlS, 

JISP, 

K, 

KAPPA, 

KSC, 

5LENGTH, 

MU  , 

MUO, 

NS3, 

CMEGA, 

P, 

PI, 

6  RMQ, 

SIGMA, 

YOR, 

YOS, 

YIP, 

YIRP, 

YIS, 

7  VISP 

REAL  JQS  ,J1S,  JISP,  JOR.Jl'^  »J1'’R,NSQ,<,>1UC 
peal  IRZl, IRZ2,  KAPPA, KSC.LENGTK, MU 
CALL  SESSIN 
o=RMC* JlSF»YlE/SIGMA 

C=  (NSQ-l.)  •K’EET  A»J  IS/  (A»PhO»SIGMA»  SIGMA) 

P=C»J1R  t  C=C»Y1R 

0  =  =’H0*  J1SF*J1R  !P  =  J1S*J1RP  t  G=  JIS'YIRP 
C='’/SIGMA 

PCG=  <G-e  )•  (G-P)  ♦  (E-0)  •  (E-D)  ♦C»CtF*F 

FCC-=FOG/  (  »G-MSQ*3)  •  (G-NSO  •?  )  *  1 E-NS C»0)  •  (E-nSC*0)  *C»C*F*pJ 
FCG  =  SORT (MUO’FOG/PPSO) 

C  ♦CWER  CALCULATICN 

DC  5  1=1,2 

FAC=SORT  (I>UO/£PSO)/FnGlI) 

P  =  (G-NSQ*6  »C*FAC)  ••2KE-NSa»0»F»FAC)**? 
PrO,((;*(5.3)»p4C).*p  ♦(c*fE-0)»FAC)*»? 
P=.t2S»FI*PI*Fl*A»A»3FTA»CM£GA*EP$0*P/RHC 
IF(P.LT.O.)  P=-F 

H(T)=EDPT(  l./P)  $  Qd)  =HtI)/FCG<I) 

5  FCG(?)=-FCG{1» 

P=l. 

return  I  end 

SU9RDUTINE  OMCAC 


COMMON  A, 

BETA  , 

BETAO, 

eiNT, 

CAPA, 

CAP?, 

20HPA(?),01DA(2) 

,DKoa(2)  , 

0“0A(?)  , 

EPSO  , 

F(?)  , 

FOG(?) 

3^UNC,  Gl?), 

gamma, 

IPZl, 

IRZ?, 

JOR, 

JOS, 

M  JIR,  JIRP, 

JlS, 

JISP, 

K, 

KAPPA, 

KSG, 

5LENGTH,  MU, 

MUO, 

NSD, 

CMFGA, 

P, 

RI, 

6  RMO,  sigma. 

YOR, 

YOS, 

YIR, 

YIPP, 

''IS, 

7  YISP 

real  IRZl, IRZ 2, kappa, Ksn, length, mu 
real  OOS ,.1S, JISP , JOR, JIP ,J1RP,NSO,K,MUO 
TE=Mi  =  SIGHA'SIGHA-nSC'RhC'PMC 
TEPM1=  TEP“1»A*J0S*Y1RP/SIGHA 
S=  SIGMA»,1  J  S=l.-?./S/S 
R=pho*A  IS=2 ./R-R*NS0*R*S 

S=S» J1S»Y1P  t  TEPM1=TFPMi»S* (N$Q»RHC»RMC/SIGMA/3IGMA-1 . )• J1S»Y0= 
S=  0H£GA»EO30*RHC*SI3"A»SIGMA 
S=  (NS0-1.)*K*K»BFTA/S 

S  =  R» (SIGNA*J  OS*Y lRtRMO»JlE»Y0  =  -?. ’JIS’YIP/ A) 

DC  5  1=1,2 

5  DHDA(I)=.E»PI»PhO*(TFPMi»c( I)  ♦S’GII)  ) 

RETURN  I  end 
SUBROUTINE  OICAC 


COMMON  A, 

beta  , 

BETAO, 

BINT, 

CAPA, 

CAP?  , 

?OHPA(?),DinA(?) 

,  CK0A(7)  , 

OMDA  (?)  , 

FPSO  , 

F(  2)  , 

FOG(  ?)  , 

TFU‘'C,  G  •?»  , 

GA M“l, 

IRZl, 

IRZ?, 

JOR, 

JOS, 

1.  JIR,  JIRP, 

JlS, 

JISP, 

<, 

kappa, 

KSG, 

SLENGTh,  mu. 

“UO, 

NSO, 

omega, 

p. 

PI, 

A  CMC,  SIGMA, 

YOR, 

YOS, 

YIR, 

YIRP, 

YIS, 

7  YISP 

peal  IRZl, ip Z?, kappa, kSO, length, mu 


B-7 


dP®*. 


tv'- 


.0^ 


NSWC  TR  81-299 


9  5aL  JOS  I'.IS,  JlS'>,JO“,J19,Jli>P,NSQ,K,‘''JO 
Te®t'‘i  =  SIG»'4»J0S»J!9»o‘-0*JlS»Ja«-7.  •JIS’JIO/A 

T£  =  M5=(NS(2-l.)*K*<»3-r4/(OM£G4'EPSD»RHO»SIG>'i*SIG‘<A)»TE!;’*5 

TEPMl  =  S1GHA*SIGMA-NSC*RH0*dHC 
TER**1=TERP  1»A»J0S»J  tP>»/SIGMA 

TE0H2=2.  /RhO/A*RHC»A*NSO»  (R  hp*  A-p  .  •phO/SIGM  A/S  1  G’IA/ A) 

TEPH2  =  T?RM3«J1S»J1P1  T£ pH3=NSO* RHC» RHC/ S IGHA/ S IGM A-1 . 
TEPflS  =  7£9M3*J1S*JCP 
TE°**1=  TEP«1  ♦TSPM?  »  TERI13 
DO  5  I=lfE 

5  DIOACI)  =-.5»Pl»RHO»(TEPHl»FJ I)  HERH5»G<I)  ) 

PETUON  $  END 
SUSfiOUTINE  OtfCAC 


CCHNON  A, 

SETA, 

BETAO, 

eiNT , 

CAPA, 

CAPS, 

2DHOA(2),CtOA(2) 

, OKOAf 2)  , 

OMOA  (2)  , 

EPSO  , 

F(2)  , 

FOG(  2) 

3FUNC,  GI2>, 

GAH-4A, 

IRZl, 

IRZ2, 

JOR, 

JOS, 

4  JIR,  JJRP, 

JIS, 

JtSP, 

y, 

KAPPA, 

KSC, 

ELENGTh,  NU, 

NUO, 

NSO, 

OMEGA, 

P, 

PI, 

6  PHO,  SIGMA, 

yor. 

y  OS, 

YIR, 

YIRP, 

YIS, 

7  YISP 
real  M(J 

real  IR21,  lPZ2»'<4PPA,'<SCtLEK'GTH 

BEAL  JOSiJlS,  JlSO,  jaB,JlR,Jlop,N$C,<,»1U0 

TEBYIiOHEG  A»mL»S1GHA»pho 

TEowi=3ETA*( JOS’YIR'SIGKA*  RhO*JiS»YOB-2.* J1S*Y1R/A) /TERMi 
YEB>*2=  (A  ‘J  OS»Y1RO  *?.•  J1S»Y1R/RhO/SIGHA/A-J1S*YOR/SIG'1A  ) 

OC  5  1=1,2 
5  Ot<PA(l)  = 

s  PI*RH0».5»<NSQ-1.)  •i<-</S:GHA»«TERMi»FtI»  ♦  T£RH?»G(I)  ) 

RETURM  t  END 
SLPOO'jTIN?  OHCAC 


CC“NCN  A, 

esTA,  ^ 

PETAC, 

BIST, 

CAPA, 

CAPS, 

PCHBA  (P)  ,  CICA  (  2) 

jOKCAC’)  , 

QNOAtP)  , 

EPSO  , 

F{2)  , 

F0G(2) 

3FUNC,  G»2), 

GA-na, 

IRZl, 

IRZ2, 

JOR, 

JOS, 

4  JIR,  JIRP, 

JIS, 

JISP, 

y, 

KAPPA, 

KSC, 

SLENGTH,  NU, 

MUO, 

MSQ, 

OMEGA, 

P, 

PI, 

6  RNQ,  SIGMA, 

YOR, 

YOS, 

TIR, 

YIRP, 

TIS, 

7  YISP 
REAL  NU 

real  J0S,J1S, J1So,J0O,J1R,J1RP,N$0,K,MU0 
real  kapp/,ksc, length 

S=  0NEGA»HL'*SIGPA»RH? 

TE°M1  =  8ET//S*  (SIGNA  •  JOS ‘JIR  ♦pmc*  J1  S*  JO  R -?  .  *  J1S»  J  1R/ A) 

TEBK  =A*J0S’J 1RP»2, •J1S»J1R/RhC/SIGMA/A 
TERN  =  TERN  -JIS’JOR/SIGNA 
00  5  1=1,2 
5  ONOA(IJ= 

=  -.5*PI*RhO»  (NS0-1.>*K»K/SIGMA»  (TERN*G(I)  »  TERMl’FC)  ) 

RETURN  J  ENO 

subroutine  error 

CC“NCN  ZZ«i.o) 

RENT  NO  6 

HRITEtG,?)  (ZZ(1)  ,1=1 ,401 
5  FCBNATdOEtZ.S) 

RETURN  t  ENO 

subroutine  GAhjcLV 


CC“N0M 

* , 

beta  , 

6ETA0, 

BINT, 

CAPA, 

CAPS, 

20H0A 12) 

,CI0A(2) 

,cyo»(2) , 

OMOA  t?)  , 

FPS5  , 

Ft2)  , 

POG(  2) 

3FUNC, 

GI2)  , 

GANha, 

IRZl, 

IPZ’, 

JOR, 

JOS, 

4  JIR, 

J1R0, 

J15, 

JISP, 

y, 

KAPPA, 

KSC, 

SLENGTH, 

NO, 

MUO, 

NSO, 

OMEGA, 

P, 

PI, 

P  RHO, 

7  YISP 

SIGMA, 

YOR, 

YOS, 

YlR, 

YIRP, 

YIS, 

CC“ncn 

/NIJNBER/ 

knax.nnax 

,IPT, IOTm 

1,AE,ALEN, 

,  SNITCH 

,z,oz 
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CC'<»*ON  sasoay 

CCHHCn  /ESY/  eaEsay 
CC«“0f!  /GSY/  Gassay 

01  “ENSIGN  G&Ri;iV(l),4ai?cAY(l),E4R  =  av(l) 
seat  J0,J1 
REat  L0A,NU,MU 
REat  LENGIf- 

“EaL  K,»(aER&,<SQ,HU(l,NSO 

DIMENSION  HF(y) 

Z  =  i3  I  DZ=L  ENGTh/ FLOAT  (NMax)  S  N**4  <  l  =  NHa;<  ♦  1 

oc  5  ir=i, Miaul 
call  itepate 

IF  (GAMMA  jLT.fl.)  PeT'.'PN 
U  BETAQ  =  kSC*  G4M“a»GA“MA 
KAOoa  =  SCRT (NSC*NS0-3FTA0) 

BETAQ  =  S(RT(££TA0) 

X=  K4PPA*a  J  J0=£FjQ(yj  tJl  =  EEJ1(») 

Tl=  NSQ*  (JQ/J  1-1  ./X)/Y 

x=GAMMa»A  t  call  MASXEL  (1  ,X .MF.OUM,  FO,  OLMI 
Fl=  SORT  (FRSO/MLQ  ) 

T?=  (-FQ/PF  (1 )  -1  ./XI  /X 

FAC  =  (KafPA'xapPA  ♦  GA"MA  »  GSMmA) •EETAC'a 

FAC  =  -K*a*KaFpa*A*xaDP4»A»Fi/FAC  •GaMMa»G4NMa 

0CA  =  FACMT1«T2)  tiOASO  =  EQA'aOA 

TI  =  K APFa«A»KAppA*a  • ( JO'JG* Jl> Jl)  -?.*Jl»Jl 

Ti  =  Ti»(fSO*  Mua»acasc/£osc) •K*3ETaG/<A?PA/KaPFA/xappa/<appa 

T?=  x»x»  (l.-FB/hF  (1  I  !♦?. 

T?=  <*S£TA0»T?/GaMM4/GaMMe/G4Hva/GAMMa  ’Jl^Jl 
T2=  rzid.YMUQ'eoaso/EPsai 

T3  =teETac»6ETAQ*NSC*><SC>/Ka=oa/K4PPa/<APP4/XaPoa 
T?  •=T3-(eETA0«e£TA0  »iCS  C » /  GA  My  A/G  AMMA/G  A  MM  A/G  AMM  1 
T3  =  T,3*  3CA»?  .•SOOT  (M'JQ/EPS  !!)»  J1*J  1 
P=.7  5»PI»(T1*T3*T'')  •SORT(FPSa/MUQ) 

IF(P.LT.0.(  P=-P 

CA°a  =  SQCT<1  ./PI  J  P  =  i. 

CAPO  =  C4FA»PCA 

garray (ID igamma  saarpav ( ir >=cap4  s  EARRAr ( rr)  =  cape 
2=?*CZ 
R  a  =  FN(Z) 

NN  =  Q  $  RETURN  J  ENTRY  GAMFIK-q 

NN  =  NN+1  S  gamma  =  GARRAX  (mnJ  J  CAPA  =  AARRaY(NN>  R  CAPS  =£  APR  A  v  ( ‘,  •.  ’ 
Y  ( N*l ) 

EET40=XSC*G4MMa»GA“Ma  SK APFa= SCOT (NSO* KSO- 2c T A B » 

BETao=  SCET(BETAO) 
rF(NN.GE.NM4Xll  NN=  0 
RETURN  J  ENO 
subroutine  iterate 


COMMON 

A  f 

beta  , 

betao, 

EINT, 

CAPA, 

cape, 

20M04 (2) 

,0I0A(2)  , 

CkOA(?» , 

OMoa  (2) , 

EPS  a  , 

F(2>  , 

FCG(?) , 

3FUNC, 

G  i2)  • 

GAM-a, 

IRTl, 

IPZ2, 

JO®, 

JQS, 

<•  JIR, 

JIRP, 

JIS, 

JtSP, 

Ky 

KAPPA, 

KSC, 

SLENGTm , 

MU, 

Mua, 

NSQ, 

CMEGA, 

P, 

PI, 

6  RMO, 

sigma, 

VOR, 

YBS, 

TIR, 

YIRB, 

YIS, 

7  YISP 
CCMNCN 

/NUMBER/ 

KMA  X,  NMAy 

,IRT, IPTMl 

,AE,ALEN 

,SMITCM, 

2,02 

RFaL  K.XSQ,  NSOjMUO.NSOPl, KAPPA 
DIMENSION  3(7QQ) 

INOEx=a 

JJM4x=Q 

KADPa=K»SQRT (  NSQ-1.) 

X=<ApoA»41  x=(3cj3  (Y>/?EJ1  (X»  ) /X 


N$noi=MSQ*l  ?  Xo.5»NS(3P1»  x 

GAyya  =  1.  i2T«F»o  (.x» /a  j  rgamma  =  .DGApya 
0N=0  (a  (GAMMA)  ?  AMT'i  =aES(  F-O) 
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?  DC  12  1=1,60 
CL  =ON 

♦  OGAKMA  «  ON=  OtA,GAt*Ma) 

JJ  =  ?»T  $3(J.j)=GAhm4  i  B(JJ‘’)  =  CS'  I  JJ1AX  =JJ“AX*2 

10  IF(aL*QN)  3,3,12 
12  CCNTIMUE 
GC  TC  6 


T  G4MMA  =  GAt^A  -  OGAHMA 

7  ESTrGAMMA  ♦  OGAI*^* A» CL /  < QL-QN) 

OE  =C<A,EST1  ?  IF  tAESCCDJ  .GT.ACIM)  GO  TC  f. 

GAmmA  =  EST  S  fi£TU9N 

8  Il!riEX=  INCEX  ♦!  %  I F{  INDEX.  GT. 70 )  GO  TO  6 
IF{QB*QLll,l,5 

1  DGAMMA  =  esT-GAMHA  J  ON  =  QS  S  GO  TO  7 

5  OGAMNA  =  GANNA  ♦OGANNA-ESTS  GANNA  =FST  SCL=OSTGO  TC  7 

6  PFINT  4,IN0EX,GAnNA,QN,DGAMNA,CL 

4  FCONAT  UNO  ,»TI-E  ITERATION  DIO  NOT  CCNVE RG E* , / ,  1 5 , 4E 15 . 9) 
PRINT  88, A 

PRINT  88, «  8(  JJl  ,  JJ=1,  JJNAX) 

83  FCPNAT(2E12.5) 


GAmna  =  -1. 
return  j  £N0 
function  C (A,GAMNAI 


CCHHON 

22, 

BETA, 

SETAO, 

SINT, 

CAPA, 

CAPC  , 

ZCHPA (2) 

,0I0A(2) 

,0X0At2)  , 

OMDAJ?)  , 

EPSO  , 

Ft2)  , 

F0G(7) 

3FUNC, 

G  (2)  , 

comma. 

IRZl, 

IPZ2, 

JOR, 

JOS, 

4  JtR, 

JIRP, 

JtS, 

JISP, 

X, 

kappa. 

KSC, 

5LEN6TH, 

MU, 

MUO, 

MSO, 

OMEGA, 

p. 

PI, 

6  RMO, 

SIGMA, 

TOR, 

T  OS, 

YIP, 

TIRP, 

TIS, 

7  >1SP 

real  jo,  J1  ,NE.:,X,XAP,t(SO,<APSO 
DIMENSION  HF  (  T)  . 

B£TSQ=KSO*GAMnA*GAM«A  $5ET  AO=SORT  (5ETSC)  S  GAMSQ=GAMM4»c.ANna 
XAPSO  =  VS'J’KSO-SETSO 
«Ap  =  SORT  (KAPSO) 

Xs  KAP»A  I  J0=8£J0<X)  J  Jt=  BEJl(X) 

TEPMl  =C  Jfl/Jl-1./X>  •A»GAMSO/l<AR 
X  =  GAMNA*4 

call  HANKEK  1  ,X,HF,CUM,F0,0UNJ 
TERM2  =-X»F0/NF«l) 

F4C=  (NSa»TERNl»T£PM?) • 1 TE°H 1*TERH2 ) - ( NSC^TERMl ♦TERMC) - ( TE RMl » IE PNJ 
21 

TE'»N1=  («NSQ»NSO-l.>  •KSQ-G&NSO)  ‘G  AMSO/xAPSC/KAPSO 
Q=P4C-TERN1 


RETURN  t  END 
FUNCTION  FN<Z) 

CC“M0M  /number/  KMAX,NNAX,IPT,IPTN1,AE,ALEN, SWITCH 
IF(SHITCMI  2,2,N 
FN  =  04  ♦Cl»EXP (02*Z)  t  RETURN 
FN  =  01  ♦C?*Z  <  RETURN 
ENTpY  FNP 
IF(SHITCH»  3,3,E 
FN  =  03»E7P(02*Z)  $  RETURN 
FN  =  02  $  RETURN 
ENTRY  INIT 
IF(SWITCH)  7,7,6 

01  =  (Z-AE)/ ( l.-EXP(SWITCH)  )  5 
04  =  Z  -01  S  P»If'T  1,02 

FCRMATt*OEXPON£NTIAL  TAPER  »(ZI  =  04  ♦  01  EXP< • , ElO .4 , ‘Z) *  > 
PRINT  8 

FCRM4T  (  =  01  INEaR  TAPEp  (wA=CUSF) 

01  =  Z  ?  02  =  <A£.Z)/ALEN 

return  1  fno 


02  =  SWITCH/ALEN  $  D3=D1»02 
S  OETURN 


A  =  A1  ♦  (A2-A1)Z»  I 


does 
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/ 
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